//mvmtest.cxx

#include <Rtypes.h>
#include "myvectormath.h"
#include <iostream>
#include <vector>
#include <cstdlib>
using namespace std;



Int_t main ( Int_t argc, const char* argv[] ){
  
  cout << "testing matrixvectorproduct..." << endl;
  cout << "test 1..." << endl;
  
  // creates diagonal matrix with real elements greater or equal zero 
  int fMatDim=3;
  vector<Double_t> A;
  vector<Int_t> JA;
  vector<Int_t> IA;
  for (int i=0; i<fMatDim; i++){
    IA.push_back(A.size());
    for (int j=0; j<fMatDim; j++){
      if(true){
	A.push_back(j+i*3);
	JA.push_back(j);
      }
    }
  }  
  IA.push_back(A.size());
  
  
  
  vector<Double_t> x(fMatDim);
  for( int i = 0; i < fMatDim; i++) {
    x.at(i)= i+2;
  }
  
  vector<Double_t> result(fMatDim);
  result=(matrixvectorproduct(&A,&JA,&IA,&x,true));
  
  cout << "A: ";
  for (auto it = A.cbegin(); it != A.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  
  cout << "JA: ";
  for (auto it = JA.cbegin(); it != JA.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  
  cout << "IA: ";
  for (auto it = IA.cbegin(); it != IA.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  
  cout << "x: ";
  for (auto it = x.cbegin(); it != x.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  
  cout << "result: ";
  for (auto it = result.cbegin(); it != result.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  
  cout << "end of test 1" << endl;
  cout << "end of testing matrixvectorproduct" << endl;
  
  
  cout << "testing skalarp..." << endl;
  Double_t skalarpresult = skalarp(&x, &result);
  cout << skalarpresult << endl;
  cout << "end of testing skalarp" << endl;
  
  cout << endl;
  
  cout << "testing vectorsub..." << endl;
  vector<Double_t> resultvectorsub=vectorsub(&result,&x);
  cout << "resultvectorsub: ";
  for (auto it = resultvectorsub.cbegin(); it != resultvectorsub.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  cout << "end of testing vectorsub" << endl;
  
  cout << "testing vectoradd..." << endl;
  vector<Double_t> resultvectoradd=vectoradd(&resultvectorsub,&x);
  cout << "resultvectoradd: ";
  for (auto it = resultvectoradd.cbegin(); it != resultvectoradd.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  cout << "end of testing vectoradd" << endl;
  
  
  
  cout << "testing VectorMultiS..." << endl;
  vector<Double_t> resultVectorMultiS=VectorMultiS(&x,skalarpresult);
  cout << "resultVectorMultiS: ";
  for (auto it = resultVectorMultiS.cbegin(); it != resultVectorMultiS.cend(); ++it)
    std::cout << ' ' << *it;
  std::cout << '\n';
  cout << "end of testing VectorMultiS" << endl;
  
  
  return 0;
  
}